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Abstract 

Many physical systems can be mapped onto solved 
or "solvable" models of magnetism. In this 
work, we have mapped the statistical mechanics 
of columnar phases of ideally helical rigid DNA - 
subject to the earlier found unusual, frustrated 
pair potential [A. A. Kornyshev, S. Leikin, J. 
Chem. Phys. 107, 3656 (1997)] - onto an ex- 
otic, unknown variant of the XY model on a fixed 
or restructurable lattice. Here the role of the 
'spin' is played by the azimuthal orientation of 
the molecules. We have solved this model us- 
ing a Hartree-Fock approximation, ground state 
calculations, and finite temperature Monte Carlo 
simulations. We have found peculiar spin order 
transitions, which may also be accompanied by 
positional restructuring, from hexagonal to rhom- 
bohedric lattices. Some of these have been exper- 
imentally observed in dense columnar aggregates. 
Note that DNA columnar phases are of great in- 
terest in biophysical research, not only because 
they are a useful in vitro tool for the study of 
DNA condensation, but also since these structures 
have been detected in living matter. Within the 
approximations made, our study provides insight 
into the statistical mechanics of these systems. 
PACS: 75.10.Hk Classical spin models, 87.15.Nn 
Properties of solutions; aggregation and crystal- 
lization of macromolecules, 64.70.-p Specific phase 
transitions, 87.14.Gg DNA,RNA 



1 Introduction 

DNA molecules in aqueous solution can be condensed into 
a variety of phases. As the density of DNA is increased, 
transitions occur from an isotropic hquid-like phase to 
a liquid crystal phase and finally to a crystalline struc- 
ture Within these phases, there exist configurations 
with different symmetries and molecular arrangements. 
For example, X-ray diffraction patterns of fibers of dif- 
ferent species of alkali metal salts of DNA reveal that 
the DNA is crystallized into several different lattice types 



12]. Likewise, liquid crystalline mesophases with differ- 
ent symmetries, including cholesteric, line hexatic, and/or 
hexagonal columnar phases, have been observed over a 
wide range of DNA concentrations P H El- These 
mesophases are relevant for several reasons - indeed they 
are seen in many systems, such as bacteria, viruses, and 
mitochondria[Hj . Understanding the structure of DNA 
aggregates may also aid in understanding the physics of 
DNA packing into sperm and phage heads and gene 
therapy vesicles [HUHIES- Last but not least, studies of 
these structures may shed light on the laws of DNA-DNA 
interaction, important, e.g., in the problem of recognition 
of homologous genes 11 . As already mentioned, the spe- 
cific phase of a DNA assembly depends heavily on DNA 
concentration, but many other factors, such as monova- 
lent salt concentration and the effects of polycationic con- 
densing agents in the solution, will also greatly infiuence 
the phase structure for a given DNA density [71^]. These 
many factors, along with the complex chemical structure 
of DNA, complicate theoretical studies of these phases 
and the transitions between them. In many studies (see, 
for example, Ref. DNA molecules are treated as uni- 
formly charged cylinders since DNA is a polyelectrolyte 
that dissociates in solution. This approximation works 
well at large interaction distances where counterions may 
screen the specific charge pattern of the DNA surface. 
But at the smaller separations where liquid crystalline 
structures are observed, a theory of the electrostatic inter- 
actions must take into account the discrete helical struc- 
ture of DNA. Recent theoretical studies JHl ^1 ^| have 
demonstrated that a number of phenomena can be ratio- 
nalized with the help of such a theory. Indeed, after disso- 
ciating in solution, DNA preserves its double helical struc- 
ture, with negative charges residing on phosphate strands 
and specifically adsorbing counter-cations settling in the 
grooves between the phosphates. This results in a heli- 
cal charge separation motif along the DNA surface which 
dictates new interaction laws at close range. Ref. [TC] 
obtained a solution for the pair potential between helical 
macromolecules in parallel alignment, having in particu- 
lar revealed that the interaction depends on the relative 
azimuthal orientation of the molecules about their long 
axes (negative charges on one molecule would like to be 
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closer to the positive charges on the other). Recent work 
has shown that this azimuthal dependence yields a rich 
phase structure for columnar aggregates and may explain 
distortions in the hexagonal columnar phase [201 112] i as 
well as a cholesteric to hexatic transition [JHl- The find- 
ings of Ref. allows a mapping of the pair interac- 
tion of DNA onto an XY-spin model of magnetism, al- 
beit with an unusual spin coupling. In this report, we 
consider the effects of temperature on the structure of 
columnar phases in order to generate their full statisti- 
cal mechanical description. Previous work [20] obtained 
ground state configurations for the columnar phases, but 
incorporating temperature effects can change the prop- 
erties of the transitions between these phases. Further- 
more, they may give rise to more phases associated with 
the relative azimuthal orientation of the DNA molecules. 
We first develop a theory of columnar assemblies fixed 
on a hexagonal lattice. We find additional Berezinskii- 
Kosterlitz-Thouless-like and topologically-related transi- 
tions in the 'spin'-structures when including temperature. 
Also, for a more complete picture of the columnar assem- 
blies, we treat positional restructuring of the molecules 
at finite temperatures. Incorporating spatial degrees of 
freedom into a Monte Carlo simulation has yielded new 
insight into anomalous spatial correlations in columnar 
phases, which may also have relevance for transitions from 
columnar to cholesteric phases. These studies reveal lat- 
tice types similar to those obtained for the ground state 
[20] but great care must be taken due to the limitations 
of the interaction potential for certain situations. For 
example, at extremely large densities where the DNA is 
closely packed, water, specifically its temperature depen- 
dent dielectric constant, can no longer be treated using its 
bulk properties as had been done to derive the interaction 
potential. Nevertheless, these results can be considered 
valid over a wide range of physical parameters and act 
as a first step to a completely atomistic approach which 
may explore a broader range of the physically relevant 
parameter space. 



properties of the solution, and they decay exponentially 
with interaxial spacing between the DNA |TB1|2S]; for the 
most updated version of these expressions, see Appendix 
A of Ref. [2n] ■ (t>i characterizes the azimuthal orientation 
of the middle of the minor groove of the i-th molecule 
relative to the direction of interaxial separation (Fig. 1). 

At large interaxial separations, the ai term dominates 
the 02 term so the ground state energy is minimized 
when the two duplexes have the same azimuthal orien- 
tations, i.e., their "spins" look in the same direction, 
(f) — (f)i ^ (f)2 = 0. But at separations below the point 
where oi — 4a2, the energy is minimized when 0^0 
and is degenerate: (p — ±(/)*, where (j)^ — arccos [ai/4a2]. 
In a hexagonal lattice, this gives rise to frustration be- 
tween the neighboring spins which results in different spin 
phases for the ground state, depending on the relative 
strength of the a coefficients. 

Since these electrostatic coefficients decay rapidly for 
increased spin separations, we need only consider nearest 
neighbor interactions. The spin-dependent term in the 
Hamiltonian for a hexagonal phase of rigid DNA frag- 
ments of length L is given by the 2D hexagonal lattice 
XY model where there is an additional frustration term 
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with the summation over only nearest neighbors. The 
spin configuration of the ground state for this Hamilto- 
nian is ferromagnetic, i.e., all the spins are aligned, if 
fli > 4a2. In the reverse case, the spins in the ground 
state are aligned in a three-state Potts [27] type configu- 
ration where the differences between the angles about any 
triangular plaquette on the lattice possess values related 
byjli 

4'2 - (t)! ^ 4>1 - 4>3 ^ 4>potts (3) 

where 



2 Hexagonal columnar assemblies 

The calculation of the pair interaction between two DNA 
in parallel juxtaposition JS] regarded DNA as consist- 
ing of a double helical charge pattern of negative charges 
along the phosphate spine of the DNA and positive 
charges adsorbed into the grooves or on the phosphate 
strands (see Fig. 1). For identical rigid helices, the 
ground state electrostatic interaction energy between two 
DNA duplexes of length L is given by 

-Eint = L[ao~ ai cos {cf)i - (/)2) + a2 cos (2 (0i - 02))] 

(1) 

where the a coefficients depend on a variety of factors such 
as the charge distribution on the DNA and the dielectric 
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Utilizing a self-consistent Hartree-Fock approximation 
(HFA) 22i ! the free energy for a given spin configuration 
can be obtained at finite temperatures. Derivations of 
these free energies are left to the appendices. For the 
ferromagnetic state where all the spins are aligned, the 
free energy is given by 
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Figure 1: The charge distribution of a single double helix (a) and the horizontal cross section of two identical double 
helices in parallel juxtaposition (b), separated by interaxial spacing R. The DNA double helix shown is considered 
to consist of two spiralling negative phosphate strands with specifically adsorbing cations in its minor and major 
grooves. The pitch H of the helix for B-DNA is approximately 34 A. In (b), (ps (which is about OAtt for B-DNA) 
is the angular half-width of the minor groove between the phosphate strands. 'Spins' characterize the aziniuthal 
orientations of the molecules; = 0i — 02 is the angle of the relative azimuthal orientation of the molecules. 



where N is the total number of spins, Chcx — —1.265 is 
a constant that depends on the geometry of the lattice, 
and J is an effective coupling that is a solution to the 
transcendental equation 
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The free energy for the Potts state is more cumbersome 
and has the form 
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where ^potts , Vi i ^'^d 772 are given in Appendix C and are 
functions of the ratio of the coupling terms, a = J2I Ji- 
This equation is closed by the additional transcendental 



equations 
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To find the transition line between the ferromagnetic and 
Potts states, the free energies (Eqs. Q and Q) are 
equated after solving the corresponding transcendental 
equations. 

To obtain a full phase diagram, we also undertook 
Monte Carlo studies of the interaction Hamiltonian, Eq. 
lO, on the hexagonal lattice. The MC simulations were 
carried out in the standard manner using the Metropolis 
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algorithm. A lattice site was chosen at random and the 
spin at this site was given a new random orientation. The 
new interaction energy between this site and its neighbors 
was calculated. If the new interaction energy was less 
than that of the original state, the new spin orientation 
was accepted. Otherwise, the new spin orientation was 
accepted with a probability according to the Boltzmann 
factor exp {—AE/ksT) where AE was the difference be- 
tween the new and old interaction energies. After equili- 
bration in the system was reached, successive MC steps 
were used to build up the canonical distribution of the 
spin configurations to obtain thermodynamic quantities 
such as specific heat, magnetic susceptibility, etc. 

In Fig. 2(a) we display the phase diagram of the sys- 
tem in terms of the relative strengths of the a coefficients 
and the temperature. Also shown on this phase diagram 
is the location of the transition between the Potts and 
ferromagnetic states obtained from the analytical expres- 
sions of the free energies. The specific heat for the ratio of 
the couplings a2/ai = 0.6 is shown in Fig. 2(b). The two 
peaks correspond to the transitions labeled as the "TSS 
transition" and the "BKT transition" (defined below) in 
the phase diagram. 

The BKT-transition corresponds to the standard Bere- 
zinskii -Kosterlitz-Thouless transition observed in many 
2D spin systems PU]. At this transition, unbinding of the 
vortices leads to an abrupt increase of the vortex density, 
defined as {l/N) ^ {vi)^ where Vi is the vorticity at lattice 

i 

site i, which in turn is defined as (l/27r) ^ (0j — the 

A 

sum of the spin differences about a triangular plaquette 
in the lattice. Figure 3 exhibits the increase in the vortex 
density at this transition. 

The peak in Fig. 2(b) found at lower temperatures 
stems from the fact that there exist two topologically dis- 
tinct ground states [2] of the Potts configuration (see 
Fig. 4(a)) akin to that seen in purely antiferromagnetic 
spin systems !32L In Fig. 4(a) the number in the center of 
each triangular plaquette corresponds to the positive he- 
licity of that triangle, defined as the sum of the clockwise 
positive change in the spin angles over 27r as the triangle 
is traversed in the clockwise direction (not to be confused 
with the vorticity where the angle difference may have 
positive and negative values). This transition is appar- 
ent in a plot of staggered, i.e., taking into account only 
downward-pointing triangles in the lattice, positive hclic- 
ity as the temperature changes (see Fig. 4(b)). Domain 
walls between these two distinct topologies are excited 
destroying the staggered helicity order at this transition, 
which we term the topological spin state (TSS) transi- 
tion. 

Comparisons can be made between the analytical 
forms, Eqs. (5)-(8), and the MC simulations to confirm 
the validity of the two methods. Specific heats of the fer- 
romagnetic and Potts states can be calculated from the 
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Figure 3: The vortex density as a function of temperature 
for the same coupling ratio as in Fig. 2(b). At the BKT 
transition, the vortex density begins to increase due to 
the unbinding of vortices in the lattice. 

free energy expressions. For the ferromagnetic state in 
the region where 4a2/ai < 1, there is quite good agree- 
ment between the simulations and the analytic forms 
(Fig. 5(a)). As expected, deviations occur at higher 
temperatures where the HFA breaks down. Agreement 
between simulation and the analytical form for the Potts 
state, however, is not as good (Fig. 5(b)). Here, the 
specific heat obtained from the MC simulations diverges 
from the analytic form at relatively low temperatures. 
This discrepancy, as well as the fact that the Hartree 
approximation can not account for the phase transition, 
can be taken as further indication of the TSS transition. 
The Hartree approximation on its own neglects all topo- 
logical excitations Therefore, this difference starts 
being significant at relatively low temperatures since do- 
main walls may be excited quite easily. Also shown in 
Fig. 5(c) is a high temperature expansion of the specific 
heat for the Kosterlitz-Thouless vortex state given by 

Ckt ^ 3 ( {Laif + {La2f \ 
ks 2 1^ (fcsT)^ J 

o( iLaif-iLa2f \ 9 ( (Laif + {La2f \ (q-. 
(fclTp ) 16 \ (k^ J' 

which conforms quite well to the simulations. 

Now that we have constructed the phase diagram for 
the 2D hexagonal lattice with spin interactions given by 
the Hamiltonian of Eq. (O for generic values of the elec- 
trostatic a coefficients, we can use their calculated values 
from Ref. I^Hl to obtain results for 'real' DNA columnar 
phases. Since these coefficients depend on the dielectric 
constant of the solution, which is temperature dependent, 
the coefficients too depend on temperature. Furthermore, 
interactions between DNA pairs are screened by charges 
in the solution, and so these coefficients are functions of 
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Figure 2: Location (a) of the various transitions in a hexagonal lattice with nearest neighbor interactions of Eq. 0. 
The TSS and BKT transitions are also observed as peaks in the specific heat (b), here shown for a coupling ratio 
iijax = 0.6. 
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Figure 4: Two topologically distinct ground states for the Potts configuration (a). At the TSS transition, domain 
walls between the two topologically distinct regions are excited so that each state becomes equally probable. This 
can be seen in a plot of the staggered positive helicity (defined in the text) as a function of temperature (b) , shown 
again for the coupling ratio of Fig. 2(b). 
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Figure 5: Comparisons of the specific heat found from analytical forms (solid lines) to that obtained from the MC 
simulations (dotted lines). Results are displayed for the ferromagnetic state (a) for a coupling ratio 02/01 =0.1 , the 
Potts state (b) for a coupling ratio 02/01 = 1.0, and finally for the vortex regime (c) for a coupling ratio 02/01 — 0.2. 



the effective inverse Debye screening length. In a poly- 
electrolyte assembly with Donnan equilibrium |34) . this 
is given by [201 

and so there is an additional temperature dependence in 
the coefficients. Here, e is the electron charge, Us is the 
salt concentration, p is the 2D DNA density in the colum- 
nar assembly, and Z \ e\ is the uncompensated DNA charge 
(the fraction of the negative phosphate backbone that is 
not compensated by readsorbed cations). After finding 
the coefficients over a range of relevant temperatures and 
densities for a given charge distribution on the DNA and 
salt concentration, we simply map these coefficients onto 
the phase diagram of Fig. 2(a) at the corresponding tem- 
perature to find where the transitions occur. 

For ambient conditions at which columnar phases are 
observed in vivo or in vitro, the BKT transition would 
unlikely be observed. The relative strength of thermal 
energy to the interaction energy, which is on the order 
of unity at the BKT transition, would engender that ex- 
tremely high, physically unviable, temperatures would be 
necessary to reach this transition. Rather, instead of go- 
ing to high temperatures, the electrostatic interactions 
could be weakened, namely by increasing the salt con- 
centration or increasing the interaxial distances between 
the DNA molecules, so that thermal energies would be 
comparable to the interaction energy. But at such low 
densities, the DNA would no longer be in a columnar ag- 
gregate. Likewise, at large salt concentrations (approxi- 
mately ten times that of physiological levels) where this 
transition could be seen, the interaction between DNA 
would be so weak that we could not assume the DNA are 
pinned to the hexagonal lattice sites. In the next section. 



we will include the effects of thermally induced positional 
restructuring in the lattice showing that such effects cer- 
tainly cannot be neglected. 

On the other hand, the TSS transition may occur at 
temperatures, DNA densities, and salt concentrations 
where columnar structures are observed in experiments. 
At this transition, the electrostatic interaction energy 
dominates the thermal energy of the DNA so that the 
DNA is essentially immobile within the assembly. In Fig. 
6 we show the location of this transition for DNA with 
different charge compensation, i.e., the fraction of the 
negative charge on the DNA phosphate backbone that 
is compensated by readsorbed cations, and also different 
salt concentration in the solution For all of these 

cases, we assumed that 30% of the cations readsorbed 
into the minor groove and 70% into the major groove 
of the DNA. To obtain this phase diagram, we took the 
length of the DNA molecules to be Lp = 50nm, which is 
the persistence length of DNA. 

We must be careful however since this transition occurs 
at relatively small interaxial spacings where the contin- 
uum electrostatic theory, which underlies the Hamilto- 
nian of Eq. |(5J), may break down. Although these spac- 
ings are quite small, they are close to the Debye screening 
length and so the linearized Poisson-Boltzmann equation 
used to derive the interaction ^Bj may still be valid at 
first approximation. However, dielectric response of wa- 
ter in the confined intracolumnar space may be different 
than in the bulk: most likely the effective dielectric con- 
stant of water will be reduced there. Its variation from 
the bulk value will vary with the density of the aggregate. 
This would slightly shift the transition to larger interaxial 
spacings, because of an increase in the electrostatic inter- 
actions between the DNA but would not eliminate this 
transition altogether. 
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Figure 6: Location of the TSS transition for different 
charge distributions on the DNA and different salt con- 
centrations in the solution. The topologically disordered 
state lies above the curves, corresponding to lower den- 
sities. The dashed line (a) corresponds to a charge com- 
pensation of 0.9 with 0.1 M salt concentration. The solid 
line (b) and dotted line (c) correspond to a charge com- 
pensation of 0.7 with 0.1 M salt solution for (b) and 1.0 
M for (c). The effects are non-trivial, as they are driven 
mainly by the variation of k (Eq.(lO)) which changes the 
balance between the electrostatic coefficients ai and 02. 

3 Positional restructuring of the 
assembly 

So far, we have assumed that the DNA molecules are 
pinned to sites on a hexagonal lattice. However, ther- 
mal motion, especially when the interactions grow weaker, 
may distort the lattice and, hence, affect the statistical 
mechanical properties of the system. Furthermore, de- 
pending on the form of the interaction, namely the in- 
terplay among the a coefficients, the hexagonal lattice 
may not be the optimal ground state of the system. Pre- 
vious experimental and theoretical works have demon- 
strated that the 2D lattice of a columnar assembly may 
be distorted hexatic @I[221 or, almost equivalently, rhom- 
bic [20]. To incorporate these effects, we must modify our 
model for the columnar phases of the previous section. 

To build up the statistical mechanical theory for colum- 
nar assemblies including positional restructuring, the 
ground state configuration of the system is first obtained 
by performing a lattice sum of the interactions among 
the DNA molecules in the assembly [20]. Upon carry- 
ing out the lattice summation, we find that the ground 
state configuration for certain values of the electrostatic 
coefficients is no longer hexagonal but rather rhombic as 
seen in Ref. [20] ■ For this structure, the spins possess a 
quasi-antiferro-magnetic-Ising (QAF) ordering: spins at 
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Figure 7: The rhombic QAF state showing the distor- 
tion from the hexagonal lattice (w > 60°). Parallel spins 
lie at opposite corners of the rhombi whereas spins with 
parameter-dependent differences lie at adjacent corners. 



opposite diagonals of the rhombic cell are the same but 
there is a finite difference between the spins at adjacent 
corners (see Fig. 7). This is quite understandable qual- 
itatively: parallel spins will repel each other due to the 
frustration-inducing a-i interaction term thus distorting 
the hexagonal lattice. The degree of this distortion again 
rests on the interplay of the a terms in the interaction. 

As before, we have derived an analytical expression for 
the free energy for this configuration taking into account 
only the six nearest neighbors. Here, however, due to the 
distortion, the electrostatic coefficients each take on two 
values according to the direction taken along the lattice: 
four of the nearest neighbors are at adjacent corners and 
the remaining two at the opposite corners along the short 
diagonal. Note, that this neglects neighbors across the 
long diagonal, but excluding these can be justified as long 
as the distortion angle is not much larger than 60° , which 
is generally the case. The free energy for this state is given 

by 
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where a = J2/J1, the form of is given in Appendix 
B, and the a coefficients are calculated at interaxial sep- 
arations of Ri, the distance between differing spins, and 
i?2 = Ri -s/S — 2cosw, the distance between parallel spins 
across the short diagonal. The equation is closed by the 
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Figure 8: The free energy for the QAF rhombic state 
(dotted hne), Potts configuration (squares) and ferromag- 
netic hexagonal states (sofid Hne) as a function of the 
average interaxial separation between nearest neighbors. 
The data are given for one persistence length long DNA 
and with charge distributions on the surface of 30% read- 
sorbed cation charge in the minor groove and 70% in the 
major groove for a charge compensation of 0.9 (below the 
break in the y-axis) and of 0.7 (above the break). For 
a charge compensation of 0.9, the QAF rhombic state is 
the minimum energy configuration below a separation of 
about 25. OA [the Potts state is the lowest energy config- 
uration at smaller spacings (not shown)], and the ferro- 
magnetic hexagonal state is the favorable configuration 
above this point. Likewise, for a DNA charge compensa- 
tion of 0.7, this transition occurs at a separation of about 
25.71. 



additional equations 
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Solving Eqs. (11) and (12), we can find the amount 
of distortion that minimizes the free energy. Figure 8 
shows comparisons between the analytic forms for the 
free energy at a temperature of 300 K for the various 




22 24 26 28 

Average interaxial separation (A) 

Figure 9: Comparisons of the level of distortion (values 
of the rhombic angle uj) obtained from the analytic forms 
(lines) and the ground state lattice summation (symbols) 
as a function of average interaxial separation. Results are 
shown for a DNA charge compensation of 0.9 (solid line 
and squares) and of 0.7 (dotted line and triangles). The 
transition from the QAF rhombic state to the ferromag- 
netic state is of first order in the HFA calculation and 
occurs at smaller interaxial spacings as compared to the 
smooth crossover obtained from ground state lattice sum 
calculations. 



lattice types for a couple of different charge distributions 
on the DNA surface. Incorporating temperature via the 
HFA calculation leads to an abrupt first order transition 
between the QAF rhombic state and the ferromagnetic 
hexagonal state, found from comparing the results of Eq. 
(5) and Eq. (11). This may well be an artifact of the HFA 
calculation [221, but nevertheless, the transition becomes 
much sharper as compared to the smoother crossover 
found from the ground state lattice sum calculation. This 
transition also occurs at a smaller interaxial spacing than 
where the ground state lattice sum calculation yields the 
crossover between the two states. Including temperature 
also alters the amount of the distortion in the rhombic 
state, as shown in a plot of w as a function of average 
interaxial spacing in Fig. 9, as compared to that found 
from the lattice sum calculations. At very large densi- 
ties, the system returns once again to a hexagonal state 
with the Potts-type spin configuration. Also, as seen in 
Fig. 8, the energy difference between the hexagonal Potts 
configuration and the QAF rhombic state is on the order 
of ksT, and so there may be a mixture of these phases 
at certain densities. Again, these results are subject to 
the limitations of the derived pair interaction at small 
interaxial spacings mentioned previously. 

As before, we perform Monte Carlo simulations on the 
2D columnar system. As opposed to the fixed lattice cal- 
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culations of the previous section, here, however, a MC 
step either corresponds to a change in the spin or a change 
in position to probe both spin and spatial degrees of free- 
dom. Again, the standard Metropohs algorithm is em- 
ployed to build up the canonical distribution of the as- 
sembly at a specified temperature. The new energy is 
found by calculating the interactions of the molecule only 
with others that lie within a specified neighborhood of the 
original molecule. This is done to increase the efficiency 
of the simulation and is justified since the interactions de- 
cay quickly with increasing intermolecular distance. From 
the simulations, we obtain relevant thermodynamic quan- 
tities, upon thermal equilibration of the system, as well as 
information concerning spin and positional correlations. 

The QAF rhombic to ferromagnetic hexagonal transi- 
tion shown in Figs. 8 and 9 also appears in these sim- 
ulations at the same average nearest neighbor interax- 
ial spacings found from our analytical explorations. In 
Fig. 10, the distribution of the spin difference between 
nearest neighbors is shown for interaxial spacings just 
below and above this transition. As shown in this figure, 
the nearest neighbor spin configuration develops from the 
two-state QAF phase to a broad single maximum about 
zero, the ferromagnetic phase. This is likewise observed 
in the spatial distribution between neighbors. Below the 
transition, the distribution shows nearest neighbor, next 
nearest neighbor, etc. correlations that match a rhom- 
bic structure, while above the transition, the distribution 
matches that of a hexagonal lattice. Note that just be- 
yond this transition at lower densities (at least for the 
charge distribution on the DNA used in the simulation 
that generated the results of Fig. 10), the DNA precipi- 
tates out of solution leaving a coexistence region of DNA 
aggregate and DNA-free solution PU] . 

In order to avoid the situation where DNA condenses 
out of solution at lower densities solely due to the choice 
of the charge distribution on the DNA, a charge distri- 
bution can be used where the interaction energy leads 
to a purely repulsive force between molecules. One such 
charge distribution is that where there the readsorbed 
cations are shared evenly between the major and minor 
grooves of the double helix. For this charge distribution, 
the interaction energy found by lattice sum calculations 
yields a flattening of the potential at intermediate den- 
sities that arises due to the spin frustration in the sys- 
tem. Due to this flattening, we flnd that as the density is 
increased, the spatial correlation function (the probabil- 
ity distribution of the location of neighboring molecules, 
A'KB?g{R) ) becomes more liquid like. Results of the MC 
simulation are shown in Fig. 11 demonstrating this effect. 
As the density is increased further, of course, the system 
once again becomes more crystalline. We essentially find 
here an effect which could be conventionally called spin- 
frustration induced melting of the positional structure of 
the columnar aggregate. 
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Figure 10: The thermally averaged probability distri- 
bution of the spin difference between nearest neighbors 
near the transition between the QAF rhombic state and 
the ferromagnetic hexagonal state. The distributions 
are shown for a DNA charge compensation of 0.9 with 
the dotted line corresponding to an interaxial spacing of 
24. SA and the solid line to the one of 2b AA at a temper- 
ature of 300K. 








Interaxial spacing (A) 

Figure 11: The probability distribution of neighbors for 
a DNA surface charge distribution with 90% charge com- 
pensation where 50% of the cations are readsorbed in the 
minor groove and 50% in the major groove. The dotted 
line corresponds to an average interaxial nearest neighbor 
spacing 24. 2A and the solid line corresponds to an aver- 
age spacing of 2Q.2A. As evident in the plot, the lower 
density distribution is more liquid-like where neighboring 
molecules have a finite probability to be at a broad range 
of interaxial spacings whereas the lower density spatial 
distribution possesses a more crystalline ordering. 
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Figure 12: Biaxial correlation functions calculated from the probability distribution of nearest neighbor spin differ- 
ences. 



In Ref. [221 analysis of the chiral electrostatic interac- 
tion for cholesteric liquid crystals demonstrated that bi- 
axial correlations, i.e., the ensemble-averaged azimuthal 
spin difference between nearest neighbors, determines the 
strength of the chiral interaction in the cholesteric phase. 
As these correlations grow larger, the chiral interaction, 
giving rise to the cholesteric phase, increases. These cor- 
relations, (cos(0)) namely and (cos(20)), can easily be 
calculated from the thermally averaged distribution of the 
spin difference between nearest neighbors in our simula- 
tion. As shown in Fig. 12, (cos((/))) obtains its maxi- 
mum value at an average interaxial separation of about 
32A at the same point that (cos(20)) becomes positive. 
Interestingly, this spacing corresponds to that at which a 
transition between the columnar phase and the cholesteric 
phase is observed in experiments |35l I36| . The transition 
between the QAF rhombic and ferromagnetic hexagonal 
phases, as shown in Fig. 9 for different charge distribu- 
tions, correspondingly occurs at this same spacing for this 
DNA charge distribution. As the interaxial spacing in- 
creases, these correlations decay to zero as a result of the 
overall weakening of the interactions between molecules 
at smaller densities, which leads to an isotropic liquid 
phase. 

4 Discussion 

We have presented results of a statistical mechanical anal- 
ysis of DNA in columnar assemblies that interact via an 
azimuthal angle ('spin') dependent pair potential |16j . 
Initially, we assumed that DNA is packed in a fixed two 
dimensional hexagonal structure so that we had to con- 
sider only 'spin-spin' interactions between neighboring 



molecules. Hence, tools similar to those used for studying 
magnetic systems could be employed. Again, this spin in- 
teraction arises from the helical charge distribution on the 
double helix. Due to a quasi-antiferromagnetic coupling 
term in the interaction, the system may be frustrated, 
which gives rise to a rich phase structure. 

Besides a Berezinskii-Kosterlitz-Thouless type transi- 
tion and an intermediate transition to a ferromagnetic- 
like state, we have found a transition associated with do- 
main formation of two distinct topologies of the spin sys- 
tem. Comparing biologically relevant values of the spin 
coupling terms to the generic phase diagram of the hexag- 
onal system, we find that this latter transition may be 
experimentally probed in dense assemblies. Calorimetric 
measurements would seem to be the simplest option to 
study this transition, but other experimental factors may 
inhibit a direct measurement of the transition in this man- 
ner. X-ray diffraction ^ or NMR techniques, how- 
ever, may be able to directly scrutinize the spin structure 
within the assembly to obtain the spin correlations in the 
system thus providing evidence if such a transition ex- 
ists. Investigations are currently being pursued by our 
experimental collaborators at Imperial. 

Extending the theory, we considered non-hexagonal lat- 
tice types and also incorporated thermally-induced posi- 
tional fluctuations of the molecules. Previous studies of 
DNA assemblies [23] demonstrated that a 2D rhombic 
(distorted hexagonal) lattice would be the ground state 
configuration of the system under certain conditions. We 
found from both analytical and numerical investigations 
a first order transition from this quasi-antiferromagnetic 
rhombic state to a hexagonal ferromagnetic state as the 
density of DNA in the system decreases. Likewise, this 
transition also appears in the probability distribution of 
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the nearest neighbor phase angle difference of the MC 
simulations. The ensemble average of nearest neighbor 
spin differences develops from a two-state spin system to 
one with a single maximum. 

Furthermore, these MC simulations revealed other in- 
teresting phenomena. For certain distributions of ad- 
sorbed cations on the DNA surface, we found that in- 
creasing the DNA density lead to a counterintuitive re- 
duction in the crystalline ordering of the system so that 
the system became more liquid like. This has been ob- 
served in classical experiments of densely packed DNA 
at roughly the same interaxial spacings as those found 
in our simulations. Again, this results from the frustra- 
tion in the spin interaction between the molecules. Also, 
in this same system, we found that the behavior of bi- 
axial correlations between neighboring molecules, which 
has been proposed as a mechanism underlying the for- 
mation of cholesteric phases at densities lower than 
those of the columnar phases, is influenced by the tran- 
sition between the quasi-antiferromagnetic rhombic and 
ferromagnetic hexagonal phases. This occurs at interax- 
ial spacings where the transition between columnar and 
cholesteric phases are observed in experiments. As the 
density of the DNA-assembly is further decreased, the in- 
termolecular interactions would grow ever weaker, until, 
finally, thermal motion would destroy any lattice order- 
ing in the system. At this point, the lattice would then 
completely melt into a liquid- like isotropic phase |l-i7l . 

Throughout this paper we have used a 2D model to de- 
scribe interactions between DNA. In three dimensions, as 
well as azimuthal angular fluctuations, there exist fluctu- 
ations of the relative positions of the molecules along their 
long axes. If the charge distribution of the DNA surface is 
taken as continuous, then these two types of fluctuations 
are indistinguishable in the roles they play in the interac- 
tion energy, i.e., an azimuthal rotation is equivalent to a 
translation along the long axis. For this case, we may sim- 
ply replace the azimuthal coordinate 4> by the coordinate 
(j) = (j) — 2ttz/ H in all our expressions, which would not 
alter our present results. However, if we consider that the 
charges on the DNA surface are discrete, this equivalence 
will be lost; there will be additional interactions which de- 
stroy this symmetry. Nevertheless, if the assemblies are 
not extremely dense, these additional interactions due to 
discreteness can be neglected and the charge pattern on 
the DNA can be considered to be continuous. Still, dis- 
creteness will restrict the corkscrew motion {(j) = 0) of 
DNA, which costs no energy for a continuous charge dis- 
tribution, within the assemblies. 

We must caution once again that the form of the DNA- 
DNA interaction tl6| we have employed is subject to cer- 
tain limitations. As the concentration of DNA in the 
solution is increased, so that the surface-to-surface sep- 
aration becomes prohibitively small, effects of nonlocal 
polarizability of the water in the narrow interstitial re- 



gions between the DNA could alter the results 38 . Fur- 
thermore, dielectric saturation and sterical constraints 
threaten to freeze the dielectric response of strongly con- 
fined water. Likewise, at such large densities, steric forces 
between DNA mediated by confined water may need to 
be also included in the interaction. Also, the pair po- 
tential was derived using a linearized Poisson-Boltzmann 
equation which is valid for weak electric fields in solu- 
tion, and thus the pair potential may not be valid when 
the electrostatic interactions are quite substantial (large 
densities). Next, as the DNA concentration is altered or 
the DNA move about, rearrangement or additional ad- 
sorption/desorption of cations on the DNA surface may 
occur in response to a change in intermolecular distances 
|84j . thus affecting the pair interaction. Last but not 
least, in this study DNA duplexes have been considered 
as ideally helical: sequence dependent distortions from an 
ideal double helical step structure of DNA j25. has been 
neglected here as well as the torsional elasticity of the du- 
plexes, that in columnar aggregates can correct the non- 
ideality EHl EH- Nevertheless, this is a first step to 
rationalize the columnar systems at finite temperatures 
before treating them with a computationally expensive 
fully atomistic approach in which these results may be 
tested. 
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Appendix A: The Hartree approxi- 
mation for the ferromagnetic state 

The partition function for our modified XY model is 

^^n/^^.'exp(-^) (Al) 

where we may rewrite Eq. (2) of the text as 

E[<P,A ^LJ: Eap(-l)ncos(p(</),,, - </),,,_!)) 

P=l 3,1 

+ C0s(p((/)j,i - + COs{p{(l)j,i - (t>j-l4 + l))] 

(A2) 

Here, we have introduced two lattice vectors Ui — jr^u 
and Vj = Irov, where u and v are unit vectors, to de- 
scribe the relative positions of the sites (DNA molecules) 
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Figure 13: Lattice labelling and relative positions of u 
and V. 



in the hexagonal lattice. The relative positions of these 
two vectors, as well as site labelling, are shown in Fig. 
13. We then start by expanding out in powers of 

(j) and dividing the energy into a Gaussian part as well as 
an anharmonic part, which contains terms of higher order 
than two in the expansion. 



£;o[0j,/] = 3NL{a2 - ai) + ^ $1 [("^^ 



+ {<P,.i - <t>j,i-i? + - (f'j-ij+i)^] (A3) 



X - 'Pj-l.lf'^ + (<t>3.l - 03-1, i)^" + (<Pj.l - 03-1. ! + 

where N {N ^ 1) is the number of sites in our system, 
and mo = L{ai — 4a2). We introduce the following lattice 
Fourier transform 



,iijku+lk^)ro 



(A4) 



where A is total area of the lattice and rp is the lattice 
spacing. We have chosen reciprocal lattice vectors corre- 
sponding to the rhombic Bravais lattice defined by ii and v 
. Then fc„ and fc„ take on values which lie within the first 
Brillouin zone for a rhombic lattice (— Tr/rg < fc„ < Tr/rp). 
This is not the only way we could choose our reciprocal 
lattice; another possible choice corresponds to the first 
Brillouin zone of the hexagonal lattice (lUI- However, we 
have been able to show the equations we obtain do not 
depend on this choice, and it is far easier to use the rhom- 
bic Brillouin zone. Using Eq. (A4) we may rewrite our 
partition function as 



Z 



D<j){k) cxp 



(A5) 



where 

Eo - ^ <^(^) (2 - 2 cos(fcaro)) 

A; 

+ {2-2 cos(fcbro)) + (2 - 2 cos((A:a - h)ro))) (l){~k) 

(A6) 



and 



ai(-l)"-^ +a2(-4)" jV 
(2n)! V" 



m = l m ■ - 



(A7) 



where km = {kmu, kmv), thus diagonalizing the 'free' part 
of the energy. 

Let us consider the free energy for Gaussian fluctua- 
tions, where we neglect Eanh- Here, we may integrate 
over 4){k) and so obtain the free energy for iV 1: 

Fo = ZNL{ai - as) 
keTNL r 



2(27r)2 



/dx / dy In ( — - [(2 — cosx) 



+ (2 - cosy) + (2 - 2cos(a:: - y))]) . 

(A8) 



The exact numerical calculation of the integral gives us 
Fo = 3NL{ai - 03) + ^^^^^ in mo + 0.235kBTNL. 
For the Gaussian correlation function, 



Jo{k) = Z-^ j D</.(fc)</.(fc)0(-fc)oxp ^- 



Eo[<t>(k)] + E^„„[<l,{k)] 



knT 



it is easy to show that 



(AlO) 



Go(k) = fcBr[mo((2-2cos(fc„ro)) 
+ (2-2cos(fe„ro)) + (2-2cos((fc.u -fc,)ro)))]~'. (All) 

To obtain the Hartree approximation we must develop 
perturbation theory for the full correlation function. We 
will omit details how one obtains the expansion j41| . 
but instead give the Feynman rules that govern the per- 
turbation expansions in our theory. On expanding out 
E[(j){k)]we need three types of diagrams (or vertices) to 
represent each term in this expansion, for n > 2 . For 
illustration some of the diagrams are shown in Fig. 14. 

Then, we may write down graphs to describe each term 
in our expansion. Each graph will contain Ny vertices, 
which are all connected to each other by lines. We assign a 
label i = 1 . . . Ny to each vertex. For each vertex i, repre- 
senting (/?", we must write down (ai(— 1)"^^ + a2(— 4)") 
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Figure 14: a.) represents all three types of vertices for the 
ip"^ term in the expansion for Eanh ■ Shown, here, are the 
four wave vectors associated with each ip'^ vertex. Also, 
we write an index r = 1,2,3 which distinguishes between 
these three types. In the corresponding expression there is 
a particular form factor associated with each of the three 
types of vertex we may draw (see text), b.) represents 
all three types of vertices for the (p^ term. In general, 
a vertex for the term is obtained by drawing a point 
with lines radiating from it. See text for a full discussion. 



^_|_ fei/(2n)! and either one of the 'form' fac- 
tors 

m=l 7n=l 

2n 

Y[ (^1 - e-'<*-"-*=""'''") (A12) 

m=l 

depending on whether the vertex is type 1, 2 or 3, re- 
spectively. Each full Feynman graph will also consist of 
Ne external lines (connected to only one vertex) each 
associated with a wave vector qj, {j = 1, . . . , Ne)- For 
each of these external lines we write down G{qj) and set 
= 'Ij O'^^ of t;he in the vertex to which the 
line is connected. There will also be Nj internal lines, 
each associated with a wave vector pfe, [k = 1, . . . , A'/), 
where each end is connected to two vertices, i and i' . For 
each of these internal lines we write down G{pk) and set 

— Kn ^ Pj foi' '^^'^ '^^ ^^'^ Kn cach of the two ver- 
tices. Then all the wave vectors for the internal lines are 
summed over. Last of all, there is also a symmetry factor 
that multiplies this, which accounts for how many ways 
a term (graph) in the expansion may be generated. 

To obtain the Hartree approximation we first consider 
only the graphs for the full correlation function shown 
in Fig. 15. The sum of these graphs we denote by Gi. 
These form a series which we may easily sum 

ksTG^^k) = Ml (2 - cos/c„ro) 

+ M2(2 - cos fc^ro) -I- M3(2 - cos(fc„ - fc„)ro) (A13) 

where 



n=0 



2mo 



+Aa2L y — 



n=0 



2kBTCr 

mo 



The zeroth order term in the series corresponds to Go . 

The first order term corresponds to the c^** vertex, the 
second order term to the (p^ vertex and so on. For we 
find the following expression: 

= N^f ? ^® " ^ cos(A:„ro))5.,i 
p 

+{2 -2cos{k^ro))5r,2 + (2 - 2cos((A;„ - kM)K3] ■ 

(A15) 

From this we may show that Cr = 1/3 and Mi = M2 = 
M3 = mi. The next step is to go back to Fig. 15 and re- 
place Go in cach loop with Gi, where Gi will be replaced 
on the l.h.s. of this expression with a new correlation 
function G2. On summing these graphs we find Eq. (A13) 
(Gi replaced by G2), but with Mi — M2 = M3 = mi 
where 



m2 



kBT\ , ( 2kBT\ 
-ai exp ( — ) -|- 4a2 exp ( — ) . (A16) 



6mi 



3m 1 



We then keep iterating this process until we have J = 
rUoo = Tioo-i and so obtain Eq. (6) of the text. 

To calculate the free energy we must consider the sum 
of graphs (Fig. 16). For the free energy in the Hartree ap- 
proximation we take care in replacing mo with J, renor- 
malizing our expansion, thus obtaining Eq. (5) of the 
text. 



Appendix B: The Hartree ap- 
proximation for the quasi- 
antiferromagnetic state 

In the QAF state we now have spontaneous symmetry 
breaking where 

{•Pj^i - (t'3-1,1) = - (l>j-i,i+i) = -0 7^ 
and(</>j,(-</>,-,;_i) =0. (Bl) 

To take account of Eq. (Bl) we rewrite — 4>j-i,i = 
ij-j + 4>'j i - 4>'j-i,i, 4>3.i - 4>]-i,i+i + - -^j-i.i+i and 
(pjd — (pjd-i ~ 't'j i~ then rewrite Eq. (A2) in 

the following form 

2 

p=i ji 

-[siii(pV)(sm(p(,^;.^, - + sin(p(</.;.,, - <?i;._i,,+i)))]] (B2) 

Expanding in powers of we may again divide the energy 
into Gaussian and anharmonic terms: 



(AM) 
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Figure 15: Graphs that are considered for the full correlation function in the first stage of the Hartree approximation. 



F = 1= - ^ 
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Figure 16: Diagrammatic expansion for the free energy. Here, Fq is the free energy calculated in the Gaussian 
approximation. 
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Eo[(f)j,i] = NL{a2 - ai) + 2NL{a2 cos(2V') - ai cos(V')) 
+ L{ai sin(V') - 2a2 sm{2ip))\(j)'f^„^^^^^y 



(B3) 



n=2 ^ -'■ 

X (0;. ( - + (ai(-l)"-i cos(V') + cos(2V')a2(-4)") 



1 

ji n=2 ^ 

where 



^4>'boundary = ^{{^'j,l-(t>'j-l,l)- 



(B4) 



mo,i = ^((ii cosV' — 4:a2Cos2tjj) and mo, 2 = -£'(01 — 4a2)- 
Again, let us consider Gaussian fluctuations and so dis- 



card all anharmonic terms. A</)^ 



boundary 



is the difl:erence in 



(p' (fluctuations in (f)) between one part of the boundary 

and another part of the boundary of our lattice, which 
cannot be neglected. We can think of the coeflicient that 
multiplies ^<P'houndary ^ ^ torque, which is non zero 
when the system is out of equilibrium. At equilibrium, 
therefore, we require that, for the ^4'boundary term to 
vanish, 

(B5) 



COS'!/' = or sini/) = 0. 
4a2 



Indeed, through our definition of ip; 



and 



must vanish, which 



occurs only if Eq. (B5) is satisfied. If only the first con- 
dition of Eq. (B5) is satisfied, only fluctuations around 
the antiferromagnetic state are under consideration. The 
second condition corresponds to fluctuations in the ferro- 
magnetic state (c.f. Appendix A). 

Now, we are able to perform the integrations over </> 
and so arrive at the free energy 

Fo = NL{ai - 02) + 2NL{ai cos(7/') - 02 cos(2V')) 
keTNL r 



+ 



2(2- 



-^-^ / dx dy In I — [mo 1 (2 — cos x) 
Try J_„ J_„ \n 



+mo,2(2 - cosy) + mo,i(2 - 2cos(a; - y))]) . (B6) 

Wc shall leave the evahiation of this integral until later. 

To obtain the Hartree approximation we proceed in 
the same manner as in Appendix A and sum the same 



graphs for both the free energy and the correlation func- 
tion. Summing up all the graphs in Fig. 15 we find that 
Eqs. (A13) and (A14) still hold with mo,i replacing mo, 
but now 



JV 



^G(p)[(2-2cos(fe„ro))5x,i 



-I- (2 -2cos{k„ro))Sr,2 + (2 - 2cos((fe„ - fe„)ro))5r,3] 

[(2 - 2cos(2;))(5^,i + (2 - 2 cos('y))(5^.2 + (2 - 2 cos(a: - y))S^,3 



(2 - 2cos(a;)) + /3o(2 -2cos{y)) + (2 - 2cos(x - y)) 



(B7) 



where /3j = mi,2/mi,i and ci = ^ C2- 

There is an important completeness relation 2ci+/3c2 = 
1. We may evaluate one of the integrals and determine 
the other through this relation. We find 



2 1 
Ci (0) = — arcsin — , 



. X 1 / 4 1 
C2(p) = 7^ 1 arcsin , = 



(B8) 



So now we are able to set Mi = M3 = mi,i and M2 = 
mi, 2 in Eq. (A13). Again, we go back to Fig. 15 and 
replace Gq in each loop with Gi, and replace Gi on the 
l.h.s. with a new correlation function G2. So we have Eq. 
(A13), but with Ml = M3 = m2,i and M2 = m2,2 where 



m2,r = —CLi exp 



kTcr{(3_ 



2m 



— ] +4a2 exp ( - 



1,0 / 



2kTCr{l3l) 



mi,o 



We again keep iterating this process until wc have 
Ji = moo,i = moo-1,1 and J2 = moo,2 = moo-i,2- Conse- 
quently, we obtain 



= 1/ |ai 



Ji = L ^ ai cos tp exp 

J2 = L 



I I ai exp ^- 



fcflT 
2Ji 
ksT 
2Ji 



ci(q) 
C2(a) ] — 4a2 exp 



4a2 cos 2i/j exp ( — — ci(a) 



Ji 



)} 



2kBT 
Ji 



C2(a) 



)} 



(BIO) 



where a = J2/J1. To obtain the free energy, again, we 
must consider the sum of graphs shown in Fig. 16. To get 
the free energy in the Hartree approximation we replace 
mo,i with Ji and mo, 2 with J2. Renormalizing the free 
energy we find 



NksT 



In 



Ji 



+ 2NL 



ao — ai cos ('(/;) exp 



keT 



2Ji 



ci(q) 



+a2 cos {2ijj) exp I — ^^^"^ ci(q) 



+ NL 
+a2 exp 



oo — ai exp 

2kBT 



Ji 

knT 
2J2 



C2(a) 



J2 



02(0;) 



(Bll) 
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where 0,af{a + ^) is determined through the following 
relationship 

/ [(2 - cos x)+a{2- cos y) + {2 - 2 cos{x - y)] \ 
^ V [(2 - cos a;) + (2 - cosy) + (2 - 2cos(x - y)] ) 

= [ da'^ I 1 - - arcsin —= ^ ] . (B12) 
Ji a'\ n ^2{a' + 1) J ^ ' 

We could not find a closed form for this integral. There- 
fore, we have approximated it with a logarithm and Pade 
approximants that interpolate between the Taylor expan- 
sion of ^af{x) in powers of x^/^ around i^a/(0) and an 
asymptotic expansion. 

Now, all the terms in -B^„^[0j,;] will contribute with 
the requirement (from Eq. (Bl)) that 

{^'o,l ~ = p™o^^ ~ exp(-ig„ro)r(g, V)) = 

{'P'j^l - <t'j-i,,+i) = jim(l - exp(-i(g„ - qy)ro)r{q,4,)) = 0. (B13) 

Eq. (B13) will constrain the value of tp as in the Gaus- 
sian approximation. We find the following expression for 

r(g,V) 

r(9, V) = 

^ ai sin(V») (-eifcj;T/(2mo,i))" - 2a2 sin(2V») (-2cifcBT/mo,i)" 

/ kBTci\ ( 2kBTci\ 

= ai sin{i/i) exp — 02 sin(2i/') exp . (B14) 

V 2mo,i / V mo,i / 

Each term in this series can be represented diagrammat- 
ically as in Fig. 17. To get the Hartree approximation, 
Eq. (B14) should be rcnormalized by replacing mo.i with 
Ji and P with a, and so replacing Go by Gqo in each of 
the loops. For Eq. (B13) to be satisfied we require that 
r(p, ip) = 0. This gives us the following equations for ^: 

ai /3fcsTci(a)\ 

cos('!/') = -; — exp or sm(t/;) = 0. (B15) 

4a2 V 2Ji / 

The set of equations Eqs. (B15) , (BIO) and (Bll) form 
a complete set for the QAF state. 

For the QAF state on a hexagonal lattice close to T = 0, 
from Eq. (B5) it follows that /3 < -1/2. The upshot of 
this is that ci(/3) and C2(/3) are complex unless /? = —1/2 
(at ai/4a2 = 1), and the QAF state at T = is unstable. 
Therefore it cannot be a ground state except at the point 
of frustration. At T 7^ we find that thermal fluctuations 
can stabilize the QAF state. However, this state always 
has a higher free energy than either the Potts or ferro- 
magnetic states, so will not be a phase of the system at 
thermal equilibrium. 

If we allow the lattice to distort to the rhombic struc- 




Figurc 18: a.) Lattice construction for the Potts state 
where and now correspond to new lattice vectors for this 
new unit cell. Here, we have constructed a new unit cell 
with the basis shown in b.). Each of the three sites in 
the basis is given a label ( 1, 2 or 3). This is to take 
into account the broken lattice symmetries of the Potts 
state. Also shown is the magnitude of the difference in 
{(p) between each neighboring site, which may either be 
ijj or 2ip. 

ture described in the text, we have 
2 

p=l jl 

+ ap{Ri) cos{p{(j>j,i - + ap{Ri) cos{p{(j>j,i - 

(B16) 

On inspection of Eq. (B16), it is easy to modify Eqs. 
(BIO), (Bll) and (B15) so as to arrive at Eqs. (11) and 
(12) of the text. 

Appendix C: The Hartree approxi- 
mation for the Potts state 

This is perhaps the most difficult of the states we must 
consider, as so many of the lattice symmetries are broken. 
To perform calculations we must make the construction 
shown in Fig. 18. Then, essentially, we must rewrite Eq. 
(A2) 

E = El + E2 + E3; 



p=i 3,1 

[cos(p(0ij,z - 4>2,j,l)) -I- COs(p(0i - 4>2,j,l-l)) 

+ COs(p(^!>ij,; - ^2,j+l,l-l)) + COs{p{^i j^i - (t)3,j,l)) 

-hcos(p((^i,j,( - 03,j,;-i)) + cos(p(^ij,; - ; 
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Figure 17: Diagrammatic expansion for T{p,ip). Here, Fq = ai sin(^/;) — 2a2sin(2^), and in the Gaussian approxi- 
mation, T{p,ip) = Fq. Usual Feynman rules apply, except there is one additional rule: the external line ending in a 
cross has no Go associated with it. 



^2 = ^EE«^'(-i)"x 

p=i jd 

[C0s(p((/)2j,i - <l>3,],l)) + C0s{p{(f>2,j,l - (/)3j-l,i + l)) 
+ C0s(p((/>2,j,i - (/-s.i-l.O) + C0s{p{(p2,j,l - 

+ cos{p{(f>2,j,i - (f>i,j-i,i+i)) + cos{p{(t>2,j,i - ^!>i,j,;+i))] ; 

[C0s(p(<^3j,; - + COs{p{4>3j^i - 

+ C0s(p((/)3j,( - (/iij+i,;)) + C0s(p((/)3j,i - (j)2,jd)) 

+ C0s(p(^3j,i - (/)2,j + l,i)) + C0s{p{<j)3j^i - (j)2,j+l,l-l))] ] 

(CI) 

where we have introduced the angles <pn where n = 1, 2, 3, 
according to their location on the lattice (the sites are de- 
fined in Fig. 18). Then wc must expand out this expres- 
sion around the Potts state, which we may do by writing 
(t>i,j,i = 4>i,j,i, 4>2,j,i = ip + 4>2,j,i and (t)3j^i =2il} + 4>i,j,i- 
The expressions we get are cumbersome and we shall not 
write them down. For brevity's sake, in what follows we 
shall only give an outline of the derivation and only state 
key results. 

The next step is to separate E into a Gaussian (£^0) 
and Enharmonic part {Eanh) in the same way as we did 
for the other states. We may Fourier transform these to 
reciprocal space using relations similar to Eq. (A4). Here 
we state a key result, that 

Eo = lY.^(-^^^o\k)m (C2) 

k 

where G(^^(fc) is given below (Eq. (C3)), = {4>i,(j)2, (ps) 
and (3i = rhi^2/'fni^i. We now define rno,i = aicos(V') — 
4a2 cos(2'(/') and TO0.2 = oi cos(2?/') — 40,2 cos{4:tJj). To ob- 
tain the Hartree result for the correlation function we 
consider the same graphs as in Fig. 15. Now, however, 
each Go is a 3 X 3 matrix, the inverse matrix of G^^ given 
in Eq. (C3). On summation of these graphs we find Eq. 
(C3), but with Gq ^, po and mo,i replaced by Gj"^, /3i and 



mi^i, respectively, where mi^i and rhi^2 are determined 
through the following relations 

"^l,s = —11 exp r +4a2 exp r 

V ™o,i J V "^0,1 J 

(C4) 

where s =1 or 2, 

'"^'^^ " 2(2^ dy [{10 + 260 + 6/3^) 

- (5 + 7/3 + 2/3^ ) (cos X + cos{y)) - (8 + 10/3 + 2/3^ ) cos(x - y) 

— /3(cos(a: — 2y) -t- cos(2a: — y))] G(a:, y) ^, 

»?2(/3) = -j^iy^ I dx I [13 + 12/3 -(4 -I- 6/3) (cos a; + cosy) 
—5 cos X cos y — 3 sin y sin x] G{x^ y)~^ (C5) 

and 

G(a;,y;/3) =6(6 + 13/3 + 6/3^) 

— 12(1 + /3)^(cosa; + cosy + cos(a; + y)) 

— 2/3(cos(a; + y) + cos(2a; — y) + cos{x — 2y)). 

(C6) 

Wc may then calculate G2^^ as wc did in the previous 
expressions. We find again Eqs. (C3) and (C4), but with 
Gj~^, mo,i, mo, 2, mi^i and mi, 2 all replaced by G^^, mi,i, 
mi, 2, m2.i and TO2.2, respectively. On further iteration (as 
described in previous appendices), we obtain Eqs. (8a) 
and (8b). Again we have used Fade approximants for 
ViiP) and miP)- 

In order to obtain the free energy again, we must con- 
sider the sum of graphs shown in Fig. 16, but with each 
Go a matrix. Then we replace mo,i with Ji and mo, 2 with 
J2, obtaining Eq. 8 with 

<i,„„,,„) = _i_,£..Y_>....(|||f}).(C7, 

r(fc, tjj) is now a matrix, and is determined by first con- 
sidering the graphs shown in Fig. 17 (where Go is also a 
matrix) and then replacing mo,i with Ji and /3 with a . 
This enables us to derive Eq. (8c) of the text. 
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ksT 



3(l + /3o) 



i(fc„-fc„)ro 



-/3o(l + e-''="''° +e-''="''») 
3(l+/3o) 



(C3) 
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